package app

// Muskingum 河道汇流：马斯京根分段连续演算法（分段马法）
type Muskingum struct {
	//========模型参数========//
	KE float64 `json:"KE"` //马斯京根法演算参数，h，敏感，KE = N * ∆t
	XE float64 `json:"XE"` //马斯京根法演算参数，敏感，0.0~0.5

	//========模型状态========//
	KL float64 `json:"KL"`
	//double KL;   //子河段的马斯京根法演算参数/h
	XL float64 `json:"XL"`
	//double XL;   //子河段的马斯京根法演算参数
	N int `json:"N"`
	//int N;    //单元河段数，即分段数
	C0 float64 `json:"C0"`
	//double C0;    //马斯京根流量演算公式I2系数
	C1 float64 `json:"C1"`
	//double C1;    //马斯京根流量演算公式I1系数
	C2 float64 `json:"C2"`
	//double C2;    //马斯京根流量演算公式O1系数
	I1 float64 `json:"I1"`
	//double I1;    //时段初的河段入流量
	I2 float64 `json:"I2"`
	//double I2;    //时段末的河段入流量
	O1 float64 `json:"O1"`
	//double O1;    //时段初的河段出流量
	O2 float64 `json:"O2"`
	//double O2;    //时段末的河段出流量
	O []float64 `json:"O3"`
	//double * O;   //各子河段出流量
	DT float64 `json:"DT"`
	//double dt;   //模型计算时段长,h
}

func (e *Muskingum) SetParams(para *Params) {
	// 默认参数
	e.DT = 24
	e.N = 1

	e.KE = para.MKE
	e.XE = para.MXE
}

func (e *Muskingum) SetState(s *State) {
	e.I1 = s.QU0
	e.I2 = s.QU
	e.O = s.O
	e.DT = s.DT
}

func (e *Muskingum) Calculate() {
	e.KL = e.DT //为了保证马斯京根法的两个线性条件，每个单元河取 KL = ∆t

	e.N = int(e.KE / e.KL) //单元河段数

	e.XL = 0.5 - float64(e.N)*(1-2*e.XE)/2 //计算单元河段XL

	denominator := 0.5*e.DT + e.KL - e.KL*e.XL //计算分母
	e.C0 = (0.5*e.DT - e.KL*e.XL) / denominator
	e.C1 = (0.5*e.DT + e.KL*e.XL) / denominator
	e.C2 = (-0.5*e.DT + e.KL - e.KL*e.XL) / denominator

	if e.O == nil {
		e.O = make([]float64, e.N) //创建存储单元流域在子河段出口断面的出流量的动态数组
		for n := 0; n < e.N; n++ {
			e.O[n] = 0.0 //单元流域在子河段出口断面的出流量为0
		}
	}
	for n := 0; n < e.N; n++ {
		e.O1 = e.O[n]                            //子河段时段初出流量
		e.O2 = e.C0*e.I2 + e.C1*e.I1 + e.C2*e.O1 //计算时段末单元流域在子河段出口断面的出流量，m3/s
		e.O[n] = e.O2                            //更新子河段时段初出流量
		e.I1 = e.O1                              //上一河段时段初出流为下一河段时段初入流
		e.I2 = e.O2                              //上一河段时段末出流为下一河段时段末入流
	}
}

func (e *Muskingum) UpdateState(s *State) {
	s.O = e.O
	s.O2 = e.O2
}
